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METHOD FOR RETRIEVING LOCAL NEAR- SURFACE MATERIAL 

INFORMATION 

FIELD OF THE INVENTION 

5 

The present invention relates to a method for retrieving 
information about near-surface material in the locality 
of a group of receivers. 

10 BACKGROUND OF THE INVENTION 

Most observations of seismic waves are made either at or 
very near to the Earth's surface. However, the elastic 
properties of the Earth close to the measurement surface 

15 show some variability. This variability is due to various 
changes in the petrophysical properties of the Earth and, 
among them, permeability changes, presence of fractures, 
presence of fluids in pores, or compaction, diagenesis or 
metarnorphism changes (Toksoz, M. , Cheng-, C- H. , and 

20 Timur, A . , 1976, Velocities o£ seismic waves in porous 
rocks: Geophysics, 41, 621-645) . It often results in data 
perturbations of at least a similar magnitude to the 
target signal ♦ 

25 Thus, before reliable subsurface information can be 
retrieved from seismic recordings at or very near to the 
earth's surface, corrections for these effects are 
required . - 
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Also, it has been shown that variations in. near-receiver 
elastic properties cause the following complications* 
Firstly, receiver static variations in the data are 
receiver- to-receiver travel time anomalies due to the 
5 propagation of most of the seismic energy through the 
heterogeneous shallow structure. Secondly, lateral 
variations in free-surface reflectivity cause both 
differences in the amount of reflected and converted 
energy, and focussing or defocusing of seismic energy. 
10 These effects result in amplitude perturbations, 
especially on horizontal recordings (KShler, s. , and 
Meissner, R., 1933, Radiation and receiver pattern of 
shear and compressional waves as a function of Poisson 
ratio: Geophys. Prosp. , 31, 421-435). 

15 

Decomposing the recorded wave fields into upgoing and 
downgoing P and S waves allows an analysis of said 
recorded wavefields without the effects of any free- 
surface interaction (Dankbaar, J\ W- M. , 1985, Separation 
20 of P- and S- waves: Geophys * Prosp. r 33, 970-986) . 

However, to perform wavefield decomposition, the free- 
surface, reflectivity and, hence, local sub-receiver 
properties, need to be known. 

25 

PRIOR ART 

In medical imaging, several concepts to estimate ^local 
material properties have been developed. In fact, medical 
30 practitioners aim to estimate local material properties 
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since quantization of mechanical properties of tissue 
could improve early detection of pathology { Ga.o, L . , 
Parker, K. J., Lerner r R. M. , and Levinson, S- F. , 1996, 
Imaging of the elastic properties of tissue - a review; 
5 Ultrasound in medicine & biology, 22, 959-977) . Unlike in 
conventional seismic surveys, the displacement field is 
measured throughout a tissue using ultrasound or magnetic 
resonance imaging (MRI ) based measurement techniques 
(Muthupillai, R. , Lomas, D, J., Rossman, P . J., 
10 Greenleaf, %J, F. , Manduca, A., and Ehman, R. L*, 1995, 
Magnetic resonance elastography jby direct visualization 
of propagating acoustic strain waves: Science, 269, 1854- 
1857) „ Quantitative elasticity reconstruction is achieved 
either by comparing modelled stress to measured strain 
15 (Gao et al w 1996, previously cited; Van Houten, E. E . 
W. , Paulsen, K. D. , Miga, M. I., Kennedy, F. E., and 
Weaver, J\ B. , 1999, An overlapping subzone technique for 
MR-based elastic property reconstruction: Magnetic 
resonance in medicine, 42, 779-786) or by direct 
20 inversion of the observed displacement field. For 
example, it has been shown that elastic properties can be 
estimated by localized inversion of the equation of 
motion (Romano, A. J., Shirron, J\ J". , and Bucaro, iT- A. , 
1998, On the non-invasive determination of material 
25 parameters from a knowledge of elastic displacements : 
theory and numerical simulation: J EES Transactions on 
Ultrasonics, Ferroelec tries , and Frequency Control, 45, 
751-759; Oliphant, T. , Mahowald, J\ _L . , Ehman, Rs, Z,. , and 
Greenleaf, J*. F. , 1999, Complex-valued ^quantitative 
30 stiffness estimation using dynamic displacement 
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Jijeasuremenfcs and local Inversion of conservation of 
momentum: IEEE Ultrasonics Symposium, pag-es 1541-154d; 
Oliphant, 2\ JEuaduca, A., E'iMan, 22. £,,, and 

Greenleaf, J. P., 2001, Complex-valued stiffness 
5 reconstruction for magnetic resonance elastography by 
algebraic Inversion of the differential equation.- Magn. 
Xeson. Med., 45, 299-310). 

Since conventional seismic data are acquired only at the 
10 surface, the problem to determine seismic subsurface 
properties is ill-posed. 

A method for estimating near-surface material properties 
from seismic data has been proposed. According to this 
15 method, the wave field is recorded using a dense 3-D 
receiver group to allow the computation of temporal and 
spatial wavefield derivatives using finite- difference 
operators. The dense 3-D receiver group permits to better 
constrain seismic near-surface velocities {Curtis, A., 
20 and Roberts son, J. o. A., 2002, Volumetric wavefield 
recording and wave equation inversion for near-surface 
material properties: Geophysics, 67, 1602-1611). The 
receiver group geometry consists of a single buried 
three-component geophone and several surface geophones. 
25 The surface geophones are sufficiently close that spatial 
wavefield derivatives can be computed, which is required 
to invert the epilation of motion for local material 
parameters-- Such geometry, which accomplishes P/S 
separation, has been originally introduced by Robertsson 
30 and Muyzert (Robertsson, J. O. A., and Muyzert, E., 1999, 
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Wave field separation using & volume distribution of three 
• component recordings: Geophys . Res. Lett., 26, 2821- 
2824) . 

5 The material properties are constrained by algebraic 
inversion of the wave equation (Oliphant et a.1 . r 1999, 
Oliphant et al f 2001 , and Curtis and Roberts son, 2002 , 
previously cited) . Denoting the recorded particle 
velocity v by [Ui u 2 u 3 ] T , the wave equation, in a 
10 homogeneous and isotropic medium, reads: 

~ = ft J V(V.v)-fl 2 Vx(Vx v), ' (1) 

where a is the P-wave velocity, p is the S-wave velocity. 
15 Boundary conditions state that the traction 013 vanishes, 
at the free surface, except for the air wave. Using the 
constitutive equation for a homogeneous elastic medium, 
this allows a substitution of first-order depth 
derivatives by expressions with horizontal derivatives: 

20 

It has been shown that the following system of equations 
is obtained after substitution of the free surface 
25 derivative conditions into the wave equation {Curtis and 
Rob&rtsson, 2002, previously cited) : 



(2a) 
(2b) 

(2c) 
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=/?M,(e)-^B,(t), 

0* _ 

dtt'o-i = frA 2 {t) - 73 

= «M 3 (*) - ;S 2 B 3 (t). 



(3a) 

(3b) 
(3c) 



The expressions for the measurable coefficients A t and B± 
5 of equations (3a) , (3b) and (3c) for direct wave equation 
inversion are obtained by substitution of the free 
surface derivative conditions,, equations (2), into the 
wave equation (1) : 



10 



M*) = ■^z(9i^ + a' 3 n) + v%v i + 2d l (v 1 f^H), (29) 

M*) = -rz (ftttt + + Vj^es 4- 25^ (V ff -vn), (30) 

4*W = ^(Vw-Vh + c^,)~V^,, (31) 

./?!(/) = 2di'(vVv / ,), (32) 

J%(«) = 2r? 2 (V„-v„) (33) 



and 



15 



^s(i) = -fr(V w -v // )-2(V? / ^) 



(34) 



The finite-difference (FD) first-order derivatives in 
depth are denoted 3jv and are given by: 



(35) 



20 
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wherein is the depth beneath the free surface of the 

buried geophcne . The vertical derivative can not be 
measured exactly at the free surface but at depth Az/2 . 
The second depth derivative 833V is given by: 



^v(A.-/4) = ' WA y >V((0 + °<^/4). (36) 



The vertical derivative at the free surface is obtained 
by substitution of the free surface conditions. If Az is 
10 small, it is possible to assume that 5 33 v (Az/4) » 3 33 v{0) - 

The following- notation is used for horizontal wavefield 
derivatives: 

[5 V* - [d x d 2 ) r (37) 

and 



= hi/. (38) 

20 

Direct wave equation inversion has been evaluated with a 
synthetic example using the model and the dense 2-D 
geophone group shown in Figure 1. Since the wavelength 
depends on the material properties, some prior knowledge 
25 has to be available to design the inter-geophone spacing. 
In the example, the P-wave velocity, a, is equal to 600 
m/s, the S-wave velocity, fS, is equal to 200 m/s and the 
density, p, is equal to 1600 kg/m 3 . A line source, 
perpendicular to the 2-D plane, located at 200 m depth, 



[0057806 ,: 15 - J an- 03 01:52 | 



8 



emits a 60 Rz Ricker wavelet- The geophones are located 
at 50 tci offset. The spacing" between the surf ac6 geophones 
is 1 meter to allow computation of spatial wavefield 
derivatives- This corresponds to, approximately, 1/6 of 
the effective wavelength. The depth beneath the free 
surface of the buried geophone is Az = 0.4 m. 

The recorded traces are shown in Figure 2 . In this 
Figure, the left graphs display recordings on the 
horizontal component and the right graphs display 
recordings on the vertical component, for both: the 
surface geophones and the buried receiver. 

Figures 3A and 3B illustrate a problem arising when 
trying to make practical use of the constraints offered 
by equations (3a) and* (3c) . Figure 3A shows the 
sensitivity of equation (3a) to variations in (3. Figure 3B 
shows the sensitivity of eqruation (3c) to variations in 
a. In both figures, the dashed and solid curves 
illustrate the left-hand side and right-hand side of 
equations (3) using the model velocities for data 
perturbed by 36dB Gaussian noise. The grey areas indicate 
the uncertainty interval for ± 10 % variations in p 
(Figure 3A) and ± 10% in a (Figure 3B) , and compare this 
to the effect of noise. It appears that derivative 
operators are relatively sensitive to effects of noise. 
Even a small amount of noise causes errors of at least a 
similar magnitude as 10%c variations in velocities. 
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Thus, if the attractiveness of the above method is that 
it can ^oa^ with the ocmr^lets v/avef ield it rs^nirss , on 
the other hand, the computation of higher order 
derivatives, which are affected by noise {Muijs, R . r 
5 Holliger, K. , and Robertsson, J. O. A , , 2002, 
Perturbation analysis of an explicit wavefx&Jd. separation 
scheme for P- and S- waves: Geophysics, 67, 1972-1982) . 
Also, errors such as mis orientation and iaislocations of 
individual geophones significantly affect estimates of 
10 spatial wavefield derivatives. 

It has been proposed to avoid the computation of the 
derivatives by substituting a plane wave into the wave 
equation (Muijs, R. , Robertsson, J\ O. A., Curtis, A., 
15 and Holliger, K. , 2003, Near -surface seismic properties 
for elastic wavefield decomposition: Estimates based on 
mul ti -component land and seabed recordings : submitted to 
Geophysics) > However, this assumes that the first break 
can be isolated in the recordings, 

20 

It has also been proposed to avoid the computation of 
spatial wavefield derivatives by transforming the wave 
equation (1) into an integral equation (J?oirtano et al., 
1998, previously cited/ Romano, . A. , Bucaro, J\ A . , Ehman, 

25 R. L. , and Shirron, J\ J., 2000, Evaluation of a material 
parameter extraction algorithm using- T4RI -based 
displacement measurements: IEEE Transactions on 
Ultrasonics, Eerroelec tries and Ereouency Control, 47, 
1575-1581) . However, the volume integrals of this method 

30 need to be evaluated and this requires interpolation of 
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the recorded particle velocity v throughout a volume Q 



line in figure 1, and the inventors of the present patent 
application did not succeed, at the priority date of said 
5 patent application, to perform such interpolation in a 
sufficient accurate manner using a buried geophone . 

In addition, it is known from the state of the art that 
it is possible to estimate the SH-velocity structure and 

10 the quality factor in a borehole by analysing an SH 
propagator matrix in the time domain. The SH propagator 
matrix is obtained from the spectral ratio of a downhole 
record over a surface record. (Trampert, J. , Cara, M. , 
and Frogneux, M, , 1993, SH propagator matrix and Q s 

15 estimates from borehole- and surfa.ee- recorded earthquake 
data.- Geophys, J\ Internat., 112, 290-299) . This approach 
is only valid for SH waves because such wave types do not 
undergo conversion at the free surface. In order to 
estimate SV and the P velocity it is necessary to 

20 consider different wave types, for example P-SV or 
surface waves. 

This method cannot be generalized to the anelastic P-SV 
case/ since the problem of constructing the propagator 

25 then becomes underdetermined, Cho and Spencer (Cho, W.H. 
and Spencer, T.W., 1992, Estimation of polarization and 
slowness in mixed wave fields: Geophysics, 805-814) 
proposed to ^use a vertical array : . of geophones to 
construct an average anelastic P-SV propagator. Assumed 

30 is that the recorded wave field can be described as a 
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superposition of plane waves. The number of geophones 
required bo construct the anelastic P-SV propagator 

depends on this number of plane waves. For example, at 
least six geophones are required to estimate the P-SV 
5 propagator for an upgoing P-wave interacting with the 

free surface, i.e. the free surface incident P-wave is 
reflected and converted at the free surface. Since this 
method constructs an averaged propagator using an array 
of geophones,. the interpretation is not clear if the 
10 medium parameters rapidly change with depth, which is 
very common close to the free surface. Therefore, this 
approach is not suitable to determine local near-surface 
material parameters. 

15 SUMMARY OF THE INVENTION" 

Considering the above, one problem that the invention is 
proposing to solve is to carry out an improved method for 
retrieving local near-surface material information which 
20 avoids the computation of wavefield derivatives so that 
it is less sensitive to measurement errors. 

The proposed solution to the above problem is, according 
to the invention, a method for retrieving local near- 

25 surface material information comprising the steps of: - 
providing a group of receivers comprising at least one 
buried receiver and at least one surface receiver 
positioned either at or very near -the Earth surface; - - 
recording a seismic wavefield; - estimating a propagator 

30 from said recorded seismic wavefield; - inverting said 
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propagator; and - retrieving said near-surface material 
information. 

The method of the invention is applicable to the whole 
5 seismogram. It is based on estimation and inversion of 
wavefield extrapolation filters. This avoids the 
computation of wavefield derivatives and, therefore, the 
method of the invention is less sensitive to measurement 
errors* This is confirmed by noise tests. 

10 

In some particular modes for carrying out the method of 
the invention, the inversion of the P-SV propagator for 
material properties is carried out in the frequency 
domain rather than the time domain; the inversion for 
15 material properties is carried out for the surface wave 
component of the seismic signal; the propagator used is 
for an anisotropic medium, specifically a transversely 
isotropic medium. 

20 DRAWINGS 

The invention will be better understood in the light of 
the following description of non-limiting and 
illustrative embodiments, given with reference to the 
25 accompanying drawings, in which: 

Figure 1 shows a model and geophone configuration 
for direct wave equation inversion according to the state 
of the art. In this figure, the dotted lines outline a 
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volume referred to previously in the present description 
of the invention? 

Figure 2 shows synthetic traces of the particle 
5 velocity, recorded by the geophone configuration shown in 
Figure 1. In this figure, traces corresponding to the 
horizontal components are on the left and traces 
corresponding to the vertical components are on the 
right; 

10 

Figures 3A and 3B illustrate the sensitivity 
analysis of wave equation inversion according to the 
state of the art; 

15 Figure 4 shows a geophone configuration and 

parameters of a half-space model according to the 
invention; 

Figure 5 shows synthetic traces of the particle 
20 velocity recorded by the surface and buried geophone 
shown in Figure 4- In this figure, traces corresponding 
to the horizontal components are on the left and traces 
corresponding to the vertical components are on the 
right/ 

25 

Figure 6 illustrates the fitting procedure for the 
data displayed in Figure 5; 

Figure 7 are cross sections intersecting the minimum 
30 of the objective function E. 
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Figure 8 illustrates the effect of stabilization 
factor c on noise -perturbed propagator filters; 

5 Figure 9 are uncertainties in estimates for a, P and 

p for different S/N ratios and c-values; 

Figure 10 shows a model and geophone configuration 
for an experiment with near-surface low velocity layer 
10 according to the invent ion ; 

Figure 11 shows synthetic traces recorded by the 
geophone configuration displayed in Figure 10; 

15 Figure 12 illustrates the effect of stabilization 

factor c in deconvolution; 

Figure 13 are cross sections intersecting the 
minimum of the objective function E; 

20 

Figure 14 are cross sections intersecting the 
minimum of the objective function E; and 

Figure 15 are uncertainties in estimates for a, £ 
25 and p for different c-values as a function of S/JsT ratio. 

It should be noted that although the illustrations 
referred to above convey 2D receiver groups anqi 
wavefields, the invention applies equally to 3D receiver 
30 groups and wave field - 
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MODE(S) FOR CARRYING OUT THE INVENTION 

The method of the invention permits to retrieve local 
5 near-surface material information. 

For the implementation of the invention, a group of 
receivers is provided. This group of receivers comprises 
at least one surface receiver placed at or very near the 
10 Earth surface, It also comprises at least one buried 
receiver placed at a depth Az under the Earth surface. 
The group of receivers may advantageously comprise a 
plurality of surface receivers and, possibly, a plurality 
of buried receivers- If the group of receivers comprises 
15 only one surface receiver and only one buried receiver, 
then the surface receiver is positioned approximately 
right above the buried geophone . If the group of 
receivers comprises a plurality of surface receivers, 
these receivers may be placed at regular distances one 
20 from each other that permit wavefield interpolation 
between the receivers in order to reconstruct the 
wavefield directly above the buried receiver. As a 
minimum, the receivers are not spaced further apart than 
half of the horizontal wavelength corresponding to the 
25 Nyquist frecpiency in the data. Also, as it will appear 
below from the description of the invention, there is no 
requirement on the depth of the buried receiver. 

The receivers of the invention are, for example, 
30 geophones. It may be single component geophones or multi- 
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component geophones, in particular, three-component 
geophonas , In the c&Be where the near surface is fluid 
saturated, geophones could be replaced by hydrophones. 

The method of the invention is based on propagator 
estimation and inversion. It is advantageously used for a 
seismic wavefield comprising P and S waves, the 
propagator being then a coupled F-SV propagator. However, 
the propagator may be calculated for the whole seismogram 
and various temporal sections of the propagator 
representing different wave types may be inverted 
separately. 

Propagator Estimation 

According to a mode for carrying out the invention, an 
elastic P-SV propagator is estimated and a scheme to 
constrain near-surface P and S velocities based on 
waveform inversion of this P-SV propagator is 
illustrated. 

In the following description, v x denotes the inline 
direction, v y the crossline direction and v z the vertical 
direction. 

Assuming that the recorded wavefield can be written as a 
superposition of plane waves, the recorded wavefield at 
depth Az can be written in the form 

v(* . J\ A.:) = P(* ? * v{*, jr, 0). (7) 
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where 

/Pn Pi2 PvA 
P - [Pn A 2 fta • (8) 

5 

Boundary conditions state that the traction vanishes at 
the free surface (z-0) and hence does not contribute to 
the downward extrapolation of the wavefield- 

For a homogeneous, isotropic and elastic medium between 
10 the free surface and depth £z, the propagator P is a 
function of the P-wave velocity a, the S-wave velocity p 
and the horizontal slowness p. For a free-surface 
incident plane wave with horizontal slowness p, the 
analytical propagator coefficients read (Afci, K., and 
15 Richards, P. G. , 2002, Quantitative seismology, 2nd. Ed.: 
University Science Books) : 



Ai (* . Pi =dV Of (f : p) + [( 1 - 20 V)/2] Gf&p) t (9) 

p T2 (t tP ) =Of(i :P ) 7 (JO) 

A»(*.P) =[(l~ Sy?V)/2] Gf (M (1 0 

Aa(*.p) = R>(I ~ 2^V)/C^.p)] Of +^MsGf(*,p), (12) 

At(*-P) - - fp9pG?[t.p) + |>(1 - 2/?V)/(2 te )] Gf (13) 

A2(«-P) =Ai(*,p) - /W-p) = Aa(i-p) = 0, (14) 



20 where 
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Gi(t, P) = [W + 9p^-) + ~ «p^s)J . (15a) 

Gf (/;, p) = [tf(* + - S(f. - qp&z)) , ( I5b) 

Gf(t, p) = [S(t + qgte) + ${t - q s Az)) , (1 5c) 

G$ {t. p) = [Sit + q s Az) - 5{t - q s As)} . ( 15d) 

These functions contain the phase shifts for two-way- 
extrapolation of the wavefield towards depth flz. The 
vertical slownesses q p and q B are given by: 

(16a) 

^=(.^ 2 -p 2 ) 1/3 . (16b) 



The coefficients of the propagator matrix [equations (9)- 
(14)] can be interpreted as follows. The amplitudes of 
the propagator coefficients are wavefield decomposition 
filters: before extrapolation of the recordings to depth 
Az, the wavefield is decomposed into upgoing and 
downgoing P, SV and SK waves. The phase shifts describe 
the two-way wavefield propagation to depth Az* Finally, 
summation of the extrapolated decomposed wavefield 
renders the total wavefield at depth Az {Qsen, A . , 
Amundsen, L., and Reitan, A. , 1999, Removal of water- 
layer multiples from multi component sea-bottom data: 
Geophysics, 64, 838-851; Aki and al . , 2002, previously 
cited) . 

The extrapolation filters for elastic P-SV wave 
propagation can be obtained directly from the data 
exploiting symmetry properties of the filters: the 
analytical expressions for the extrapolation filters show 
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that P llr P 22 and P 33 are even function around t = 0 and. Pi 3 
and. p 31 are odd functions. Hence, the spectra cf P llf p 2 2 
and P33 are entirely real, whereas the spectra of Pi 3 and 
P31 are purely imaginary. These properties are used to 
estimate the extrapolation filters without prior 

information on oc, p and p. E<juatingr real and imaginary 

parts of equation (7) in the frequency domain shows that 

the propagator coefficients are given by: 



10 



In these equations, 9* [ Wa>, sr) I denotes the real part of 

\t{<d,x,z) and 3[v((d,x,z)] is the imaginary part of 

vi<x>,x,z)* Note that these data-estimated propagator 
15 filters do not require an' assumption that only one single 
slowness event exists in the data. Expressions for P(x,t) 
are found £?y calculating the inverse Fourier transform of 
epilations (17) to (21) ♦ 

20 The propagator terms P U/ P 22 and P 33 have the greatest 
amplitudes and, especially, for waves with near-vertical 
incidence, the off-diagonal terms P 13 and P 31 of the 
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propagator are insignificant as compared to P u , P22 and 
P 33 which are mainly used to constrain the P and S waves 
velocities. Then, approximations are possible to obtain 
Pn and Pn can be approximated using only the 

5 horizontal components v K , both at the surface and at 
depth. P33 can be approximated using only the vertical 
components. As a result, a group of receivers comprising 
only one single- component surface geophone and only one 
single- component buried geophone placed right above the 
10 surface geophone may be used for the implementation of 
the invention r The phase information is not changed by 
the approximations made. However, the amplitude 
information is distorted to some extent. 

15 Spectral divisions in equations (17) to (21) are 
numerically unstable because signals are band-limited and 
due to the existence of notches in the spectrum. The 
following procedure is used to stabilize said spectral 
divisions. For arbitrary signals F(«) and S(gd), the 

20 spectral division of F(cd) by S (<o) is given by: 

In practice, the spectral ratio is estimated using the 
25 following technique (Langston r C. A., 1979, Structure 
under Mount Rainier, Washington, inferred from 
teleseismic hody waves; J\ Geophys. Res., 84, 4749-4762): 
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where 

5 *„ = iiinx{5(£«/)5 + (u;),cuiax [5(a;)5*(w)]} (24) 

and W(ca) is a frequency window to limit the final 
frequency band in the estimated deconvolution , The 
complex conjugate of S is denoted by S* . The function $> 33 
10 can be thought of as simply being the autocorrelation of 
S(co) with any spectral troughs filled to a level 
depending on the parameter c . The frequency windowing 
function W(co) is a tapered window with cut-off 
frequencies the minimum and maximum frequency for which 

15 

. S(*)£*M > c max [S(*>)$*{ur)] . (25) 

This criterion implies that the parameter c also controls 
the bandwidth of the spectral ratio G' (tf) . 

20 

Propagat or Inver s i on 

The inversion scheme may be illustrated for near-surface 
velocities with a half -space example. In Figure 4, The P 
25 velocity is 600 m/s, the S velocity is 200 m/s, and the 
density is 1600 kg/rn 3 . An explosive line source is located 
at 200 m depth and emits a 120 Hz Ricker wavelet. Multi- 
component geophones are positioned at 50 m offset, one is 
located at the free surface and the second geophone is 
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located at 1,0 m depth. Synthetic data are computed using 
the Cagniard-de Hoop method (de Hoop, A. T. r and van der 
Hljden, J\ H. M. T. , 1983, Generation of acoustic waves 
by an impulsive line source in a fluid/solid 
5 configuration with a plans boundary: J\ Acoust. Soc. Am., 
74 r 333-342) , Figure 5 shows synthetic traces of the 
particle velocity recorded by the surface and buried 
geophones shown in Figure 4. The left graph shows traces 
of v x and the right graph shows traces of v z . 

10 

Both the propagators Pn, P13, P31 and P 33 calculated from 
the data displayed in Figure 5 and the fitting procedure 
to estimate cx, p> and p are outlined in Figure 6 , 
Analytical solutions for propagator filters are functions 

15 of a, p and p. Given estimates for oc, £ and p, such 
analytical solutions can be computed as shown on the 
right graph (a) of this Figure . In practice, data is 
band- limited, and these analytical solutions can not be 
compared directly to the estimated filters. So, the 

20 frequency band of these analytical expressions is limited 
to a similar frequency band as the data-estimated 
propagator filters . Graph (b) shows the time-domain 
expressions of frequency filters W(<a) [equation (23)], 
which limit the solution to the frequency band in which 

25 the spectral division was performed. A convolution of the 
analytical solution in graph (a) with the filter shown in 
giraph (b) allows a comparison to the data-estimated 
propagator. Graph (c) illustrates the fit between the 
data- estimated and band-limited analytical propagator 
30 filters. It shows that there is a good agreement between 
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the analytical (thin lines) and data-estimated propagator 
filters (thic3c 1 inss ) . The s tsi) n_ 1 x z at 2. on f set or c — 0.05. 

Figure 6 shows that both the phase and the full waveform 
5 of the extrapolation filters can be used to constrain 
near-surface properties. For example, the time delay 
between the peaks in P u and P33 gives the two-way vertical 
traveltime for SV and P waves , respectively. The 
waveforms contain both phase and amplitude information, 
10 The amplitudes of the propagator filters are controlled 
by the free-surface reflectivity. To illustrate the 
constraints of propagator waveform fitting on a, £ and p, 
the following objective function is defined: 

[5 E = E n + + i? i3 + Esu (26) 

with 



20 

The analytical solution for the propagator filter 

component ij is denoted by ^ ,a . , ^P) < the data-estimated 
propagator filter component ij is denoted Pij(t). Note 

that different components within the sum on the right 
25 hand side of equation (2 6) could be weighted differently 
to account for relative uncertainty or arty other effect, 
if desired. 



(0057806 15 Jan 03::p r 52, :| 



24 



Cross sections intersecting the minimum of the objective 
function E are shown in Figure 7. Contours are drawn for 
E = 0,10 and E - 0.20. These show that perturbations in p 
have a relatively small influence on estimates of ct and 

Figure 8 illustrates the effect of stabilization factor c 
on noise-perturbed data-estimated propagator filters for 
recorded data perturbed by 20 dB Gaussian noise. Data- 
estimated propagator filters are compared to the 
analytical filters which are computed for different 
values of the stabilization factor c [equation (24)], 
namely c = 0,02 (Figure 8, left graph), c = 0.05 (Figure 
8, centre graph) and c = 0.20 (right graph). The thin 
lines are the analytical solutions and the thick lines 
the data-estimated propagator filters. The amplitudes of 
each filter are normalized to the solution for c » 0.02. 
For increasing c-value, the effect of noise is reduced in 
the propagator filters. Note the differences between the 
different coefficients: for a low c-value, there is an 
excellent fit between and P 13 , whereas noise 

significantly affects P 31 . This is related to the absence 
of free-surface incident S-waves in this example. A 
second effect is that the frequency band is limited for a 
higher c-value. The filters are smoothed and contain less 
energy. Some care must: be taken selecting the 
stabilization factor c, since the analytical solutions 
t ep:e not corrected for the stabilization of notches in the 
spectral divisions r which have to be evaluated to obtain 
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Figure 9 shows the uncertainties in estimates of a, 6 and 
p for different S/N ratios and c-values. For example, the 
relative standard deviation a in a is given by: 



Cr(or - ftp) _ J_ 



,V- I. 



1/2 



(28) 



where oc 0 is the model P-wave velocity and «i is the 
estimated P-wave velocities. For each noise level, 
10 experiments were conducted 1000 times <N = 1000) with 
different manifestations of Gaussian noise. The minimum 
of the objective function E was determined using a 
forward search method. It appears that the estimates of 
oc, p and p are most robust for c = 0.20, and that a. and p 
15 are better constrained than p. Errors in estimates of p 
do not influence estimates of a. and p significantly. 



In the previous example only a single slowness arrival 
was recorded. The second experiment is performed in a 

20 model with a near-surface low velocity layer. 
Reverberations in the near-surface low velocity layer 
result in multiple arrivals- It is important to assess 
the consequences of multiple arrivals with different 
slownesses, since a single slowness assumption is 

25 required to construct the analytical propagator filters. 

Figure 10 shows a model and geophone configuration for 
experiment with near-surface low velocity layer. The 
receiver group is similar to the previous experiment. The 
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near-surf ace layer is 5 m thick with a = 600 m/s, (3 - 200 
m/s and p = 1600 kg/m 3 . The parameters of the underlying 
half-space are: a = 1500 m/s, p = 400 m/s and p - 1800 
kg/in 3 . The P-source is located at 100 m depth and emits a 
5 120 Hz Ricker wavelet. 

The recorded synthetic traces are shown in Figure 11* 
Recordings of v x are shown in the left graph and 
recordings of v z are shown in. the right graph. 

Figure 12 compares data estimated propagator expressions 
to analytical propagators. As for Figure 8, data- 
estimated propagator filters are compared to the 
analytical filters which are computed for different 

15 values of the stabilization factor c: c = 0.02 (Figure 
12, left graph), c- = 0.05 (Figure 12> centre graph) and c 
= 0.20 (right graph). The thin lines are the analytical 
solutions and the thick solid lines are the filters 
estimated from the data- Traces are normalized with 

20 respect to the c = 0.02 curves. It appears that r again, 
both the temporal resolution and the noise level are 
reduced with increasing c-value. However, the fit between 

and P13, for example, deteriorates for c > 0.05. This 
is caused by corrections for notches in the spectral 

25 divisions, which need to be evaluated to obtain Pij. The 
analytical solutions are not corrected for the 
stabilization of notches in the spectral divisions. Also, 
for high > c-values , the propagators lack clear maxima 
introducing several minima in the waveform fitting misfit 

30 function E [equation (26) ] . 
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Cross sections intersecting the minimum of these misfit 
functions E are shown in Figures 13 and 14 for c . = 0.05 
arid, c = 0*20. Contours are drawn for E = G _ 1 5 , E - 0.20 
5 and c - 0.05. It appears that, for c = 0.05, the minimum 
is located close to the model velocities, whereas for c = 
0.20, a low value for (5 is preferred by waveform 
inversion. 

10 Figure 15 shows uncertainty in estimates for oc and (3 for 
different c-values as a function of S/N ratio and 
illustrates the effect of noise on velocity estimation. 
It appears that propagator waveform inveraion for a and (3 
is robust if the c-value is properly selected. For low c- 
15 values, for example c ~ 0.02, noise affects the velocity 
estimates whereas, for high c-values, errors are 
introduced due to notches in the spectrum and non-unicjue 
solutions exist. Multiple arrivals, i.e. the single 
slowness assumption for propagator waveform fitting, 
20 cause only small errors in the velocity estimates. 

From the above description, it results that the method of 
the invention permits to retrieve local near-surface 
material information without imposing any constraints on 
25 the depth of the buried geophone. It avoids the 
computation of spatial wavefield derivatives and, 
therefore, it is less sensitive to noise. 

The propagator filters are estimated directly from the 
30 data. In order to retrieve the near-surf ace velocities, a 
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single slowness assumption has been made. However, in a 
medium with a near-surface low-velocity layer. the 
results are not unduly influenced by multiple arrivals, 
The synthetic examples show that the inverted P and S 
wave velocities are not sensitive to variations in the 
horizontal slowness. Propagator inversion provides robust 
constraints on local near-surface properties. 

While the invention has been described in conjunction 
with the exemplary embodiments described above, many 
equivalent modifications and variations will be apparent 
to those skilled in the art when given this disclosure. 
Accordingly, the exemplary embodiments of the invention 
set forth above are considered to be illustrative and not 
limitating. Various changes to the described embodiments 
may be made without departing from the spirit and scope 
of the invention. 

Specifically, the approach could be used on surface wave 
data rather than P-SV plane waves with a corresponding 
change in the analytic propagator used in the inversion 
Stage. 

The approach could also be used in the transversely 
isotropic case, again with an appropriate change in the 
analytic propagator ♦ 
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1. A method for retrieving local near-surface material 
information comprising; the steps of: 

5 providing a group of receivers comprising at least 

one buried receiver and at least one surface receiver 
positioned either at or very near the Earth surface; 
recording a seismic wavefield; 

estimating a propagator from said recorded seismic 
10 wavefield; 

inverting said propagator; and 

retrieving said near-surface material information. 

2. The method of claim 1, wherein the group of 
15 receivers comprises a plurality o£ surface receivers , 

3. The method of claim 1 or claim 2 wherein the 
receivers are geophones . 

20 4. The method of any one of claims 1, 2 or 3 , wherein 
the buried receiver is a three-component geophone. 

5. The method of any one of the previous claims, 
wherein the seismic wavefield comprises P and S waves. 

25 

5, The method of any one of the previous claims, 
wherein the propagator is calculated for the whole 
seismogram. 
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7. The method of any one of the previous claims, 
wherein the propagator is a coupled P-SV propagator. 

8. The method of any one of the previous claims, 
further comprising the step of: 

assuming that the recorded seismic wavef ield can be 
written as a superposition of plane waves. 



10 



15 



20 



9- The method of any one of the previous claims, 
further comprising the steps of: 

defining propagator coefficients, which are 
wavef ield decomposition filters, for the free-surface 
plane wave; and 

extrapolating said coefficients to depth flz. 

10.. The method of any one of the previous claims, 
wherein the propagator P(x, t) is obtained by calculating 
the inverse Fourier transform of the following 
coefficients : 
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where , 1R [ v{o>, z) ] denotes the real part of \r{ca, x, z) and 

3 [ v(gd, jc, z) I is the imaginary part of v{u,x r z) f v x is the 

inline velocity component and v x the vertical velocity 
component . 

5 

11- The method of any of the previous claims, wherein 
the inversion of the P-SV propagator for material 
properties is carried out in the frequency domain . 

10 12. The method o£ any of the previous claims, wherein 
the inversion for material properties is carried out for 
the surface wave component of the seismic signal. 



15 



13. The method of any of the previous claims, wherein 
the propagator used is for an anisotropic medium 
specifically a transversely isotropic medium. 
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ABSTRACT 

The invention concerns a method for retrieving local 
near-surface material information. According to the 
invention, the method comprises the steps of: - providing 
a group of receivers comprising at least one buried 
receiver and at least one surface receiver positioned 
either at or very near the Earth surface; - recording a 
seismic wavefield; - estimating a propagator from said 
recorded seismic wavefield; - inverting said propagator; 
and - retrieving said near-surface material information. 
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Fig. 10 
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